#pragma once

#include "gtest/gtest.h"

#include "../../Factory/ElementCalculatorFactory.h"
#include "../ShearElementCalculator.h"

using SG::Algebra::Matrixd;
using SG::DataStructure::Common::Real;
using SG::DataStructure::Post::SHEARForceData;
using SG::DataStructure::Post::SHEARStressData;
using namespace SG::DataStructure::FEM;
using namespace SG::FEMSolver::Element;

// ElementCalculatorFactory 是否正确创建 Shear 单元计算类
TEST (ElementCalculatorFactoryTest, ShearElementCalculator)
{
    ElementCalculatorFactory::Register (ElementTypes::ELETYPE_Shear, ShearElementCalculator::Create);

    // Test data
    auto calculatorPtr      = ElementCalculatorFactory::Create (ElementTypes::ELETYPE_Shear);
    auto shearCalculatorPtr = dynamic_cast<ShearElementCalculator*> (calculatorPtr.get ());

    // Check
    ASSERT_TRUE (nullptr != shearCalculatorPtr);
}

class ShearElementCalculatorTest : public ::testing::Test
{
  public:
    virtual void SetUp () override
    {
        // 指定材料属性
        E  = 100000;
        nu = 0.3;
        G  = 38461.5384615385;

        points        = SG::Algebra::PointArray (4);
        points[0].m_x = 0.0;
        points[0].m_y = 0.0;
        points[0].m_z = 0.0;
        points[1].m_x = 0.627461386417860;
        points[1].m_y = 0.0110581128284480;
        points[1].m_z = 0.0;
        points[2].m_x = 1.0;
        points[2].m_y = 0.982376424959660;
        points[2].m_z = 0.0;
        points[3].m_x = 0.0;
        points[3].m_y = 1.0;
        points[3].m_z = 0.0;

        // 截面属性
        property.m_thickness = 0.01;
        property.m_F[0]      = 0.5;
        property.m_F[1]      = 0.5;

        // 指定材料属性
        E2  = 210000.0;
        nu2 = 0.3;
        G2  = 80769.234;

        // 截面属性
        property2.m_thickness = 0.35;
        property2.m_F[0]      = 0.0;
        property2.m_F[1]      = 0.0;
    }

  protected:
    Real                    E;
    Real                    nu;
    Real                    G;
    SG::Algebra::PointArray points;  ///< 单元节点坐标
    PSHEARData              property;

    // used for stiffness overall test
    Real E2;
    Real nu2;
    Real G2;

    PSHEARData property2;
};

TEST_F (ShearElementCalculatorTest, computeShearCoordTrans)
{
    ShearElementCalculator  elementCalc;
    Real                    area;
    Real                    warp;
    Matrixd                 Tel (3, 3);
    Matrixd                 Tel23 (2, 3);
    SG::Algebra::PointArray localPoints (4);
    elementCalc.ComputeCoordTrans (points, Tel, localPoints, area, warp);

    Matrixd TelTarget (3, 3, { 0.999844740966936, 0.0176208388211154, 0, -0.0176208388211154, 0.999844740966936, 0, 0, 0, 1 });
    // Matrixd Tel23Target(2, 3, {0.999844740966936, 0.0176208388211154, 0,
    //                         -0.0176208388211154, 0.999844740966936, 0,});
    ASSERT_TRUE (Tel == TelTarget);
    // ASSERT_TRUE(Tel23 == Tel23Target);

    SG::Algebra::PointArray localPointsTarget (4);
    localPointsTarget[0].m_x = 0;
    localPointsTarget[0].m_y = 0;
    localPointsTarget[0].m_z = 0;
    localPointsTarget[1].m_x = 0.627558820593536;
    localPointsTarget[1].m_y = -1.73472347597681e-18;
    localPointsTarget[1].m_z = 0;
    localPointsTarget[2].m_x = 1.017155038;
    localPointsTarget[2].m_y = 0.964603063324701;
    localPointsTarget[2].m_z = 0;
    localPointsTarget[3].m_x = 0.0176208388211154;
    localPointsTarget[3].m_y = 0.999844740966936;
    localPointsTarget[3].m_z = 0;
    ASSERT_TRUE (localPoints == localPointsTarget);
}

TEST_F (ShearElementCalculatorTest, ComputeStiffness)
{
    ShearElementCalculator elementCalc;
    Real                   area;
    Real                   warp;
    Matrixd                Tel (3, 3);

    SG::Algebra::PointArray localPoints (4);
    elementCalc.ComputeCoordTrans (points, Tel, localPoints, area, warp);

    Real pSet[4];
    auto adaptability = elementCalc.computeAdaptability (localPoints, area, G, property, nu, pSet);

    Real adaptabilityTarget = 0.000437701183762857;
    ASSERT_NEAR (adaptability, adaptabilityTarget, SG::DataStructure::Common::TOLERANCE);

    Matrixd Kel (12, 12);
    Matrixd Tel23 (2, 3);
    Tel23 (0, 0) = Tel (0, 0);
    Tel23 (0, 1) = Tel (0, 1);
    Tel23 (0, 2) = Tel (0, 2);
    Tel23 (1, 0) = Tel (1, 0);
    Tel23 (1, 1) = Tel (1, 1);
    Tel23 (1, 2) = Tel (1, 2);
    Matrixd T (8, 12);
    T.BlockFill (0, 0, Tel23);
    T.BlockFill (2, 3, Tel23);
    T.BlockFill (4, 6, Tel23);
    T.BlockFill (6, 9, Tel23);
    elementCalc.ComputeStiffness (localPoints, adaptability, area, E, property, Kel);
    Kel = T.Transpose () * Kel * T;

    Matrixd KelTarget(12, 12, {
        505.208671136838,  117.376762781801, 0, -323.170367532689, -116.621871573759,  0,   -112.43623676682, -110.454708310906, 0, -69.6020668373289,  109.699817102864, 0,
        117.376762781801,  305.333705620365, 0,  61.4533751185632, -107.888505490641,  0,  -110.454708310906,  -108.50810147043, 0, -68.3754295894584, -88.9370986592938, 0,
                       0,                 0, 0,                 0,                 0,  0,                  0,                 0, 0,                 0,                 0, 0,
       -323.170367532689,  61.4533751185632, 0,  460.106112713062,  2.23433223173047,  0,  -93.8495672066756, -131.595832182603, 0, -43.0861779736976,  67.9081248323092, 0,
       -116.621871573759, -107.888505490641, 0,  2.23433223173047,  271.986268068417,  0,   46.4794145097192, -57.0677674583285, 0,  67.9081248323092, -107.029995119447, 0,
                       0,                 0, 0,                 0,                 0,  0,                  0,                 0, 0,                 0,                 0, 0,
        -112.43623676682, -110.454708310906, 0, -93.8495672066756,  46.4794145097192,  0,   383.133273352696,  169.331789008882, 0, -176.847469379201, -105.356495207695, 0,
       -110.454708310906,  -108.50810147043, 0, -131.595832182603, -57.0677674583285,  0,   169.331789008882,  273.418927932342, 0,  72.7187514846274, -107.843059003584, 0,
                       0,                 0, 0,                 0,                 0,  0,                  0,                 0, 0,                 0,                 0, 0,
       -69.6020668373289, -68.3754295894584, 0, -43.0861779736976,  67.9081248323093,  0,  -176.847469379201,  72.7187514846274, 0,  289.535714190227, -72.2514467274782, 0,
        109.699817102864, -88.9370986592939, 0,  67.9081248323092, -107.029995119447,  0,  -105.356495207695, -107.843059003584, 0, -72.2514467274782,  303.810152782325, 0,
                       0,                 0, 0,                 0,                 0,  0,                  0,                 0, 0,                 0,                 0, 0
       });

    ASSERT_TRUE (Kel == KelTarget);
}

TEST_F (ShearElementCalculatorTest, ComputeStiffnessParallCase)
{
    auto points   = SG::Algebra::PointArray (4);
    points[0].m_x = 1.0;
    points[0].m_y = 0.0;
    points[0].m_z = 0.0;
    points[1].m_x = 1.0;
    points[1].m_y = 1;
    points[1].m_z = 0.0;
    points[2].m_x = 0.0;
    points[2].m_y = 1.0;
    points[2].m_z = 0.0;
    points[3].m_x = 0.0;
    points[3].m_y = 0.0;
    points[3].m_z = 0.0;
    ShearElementCalculator elementCalc;
    Real                   area;
    Real                   warp;
    Matrixd                Tel (3, 3);

    SG::Algebra::PointArray localPoints1 (4);
    elementCalc.ComputeCoordTrans (points, Tel, localPoints1, area, warp);

    Real pSet[4];
    auto adaptability1 = elementCalc.computeAdaptability (localPoints1, area, G2, property2, nu2, pSet);

    Real adaptabilityTarget1 = 1.7687074122449008e-05;
    ASSERT_NEAR (adaptability1, adaptabilityTarget1, SG::DataStructure::Common::TOLERANCE);

    Matrixd Kel1 (12, 12);
    Matrixd Tel23 (2, 3);
    Tel23 (0, 0) = Tel (0, 0);
    Tel23 (0, 1) = Tel (0, 1);
    Tel23 (0, 2) = Tel (0, 2);
    Tel23 (1, 0) = Tel (1, 0);
    Tel23 (1, 1) = Tel (1, 1);
    Tel23 (1, 2) = Tel (1, 2);
    Matrixd T (8, 12);
    T.BlockFill (0, 0, Tel23);
    T.BlockFill (2, 3, Tel23);
    T.BlockFill (4, 6, Tel23);
    T.BlockFill (6, 9, Tel23);
    elementCalc.ComputeStiffness (localPoints1, adaptability1, area, E2, property2, Kel1);
    Kel1 = T.Transpose () * Kel1 * T;

    Matrixd Kel1Target(12, 12, {
        7067.30797000000	,-7067.30797000000	,0       	,-7067.30797000000	,-7067.30797000000	,0	     	,-7067.30797000000	,7067.30797000000	,0	      	,7067.30797000000	,7067.30797000000	,0,
        -7067.30797000000	,7067.30797000000	,0       	,7067.30797000000	,7067.30797000000	,0	     	,7067.30797000000	,-7067.30797000000	,0	      	,-7067.30797000000	,-7067.30797000000	,0,
        0	                ,0	                ,0       	,0	                ,0	                ,0	     	,0	                ,0                 	,0	      	,0	                ,0               	,0,
        -7067.30797000000	,7067.30797000000	,0       	,7067.30797000000	,7067.30797000000	,0	     	,7067.30797000000	,-7067.30797000000	,0	      	,-7067.30797000000	,-7067.30797000000	,0,
        -7067.30797000000	,7067.30797000000	,0       	,7067.30797000000	,7067.30797000000	,0	     	,7067.30797000000	,-7067.30797000000	,0	      	,-7067.30797000000	,-7067.30797000000	,0,
        0	                ,0	                ,0        	,0	                ,0                 	,0	     	,0	                ,0                  ,0	      	,0	                ,0                  ,0,
        -7067.30797000000	,7067.30797000000	,0       	,7067.30797000000	,7067.30797000000	,0	     	,7067.30797000000	,-7067.30797000000	,0	      	,-7067.30797000000	,-7067.30797000000	,0,
        7067.30797000000	,-7067.30797000000	,0       	,-7067.30797000000	,-7067.30797000000	,0	     	,-7067.30797000000	,7067.30797000000	,0	      	,7067.30797000000	,7067.30797000000	,0,
        0	                ,0                  ,0       	,0	                ,0                  ,0	     	,0                 	,0                  ,0	      	,0	                ,0                 	,0,
        7067.30797000000	,-7067.30797000000	,0       	,-7067.30797000000	,-7067.30797000000	,0	     	,-7067.30797000000	,7067.30797000000	,0	      	,7067.30797000000	,7067.30797000000	,0,
        7067.30797000000	,-7067.30797000000	,0       	,-7067.30797000000	,-7067.30797000000	,0	     	,-7067.30797000000	,7067.30797000000	,0	      	,7067.30797000000	,7067.30797000000	,0,
        0	                ,0                  ,0       	,0	                ,0                  ,0	     	,0                  ,0                	,0	      	,0                 	,0                  ,0
    });
    ASSERT_TRUE (Kel1 == Kel1Target);
}

// 1,4节点边平行 2，3节点平行边 （测试tol数值为 0.01）
TEST_F (ShearElementCalculatorTest, ComputeStiffnessS14ParallS23)
{
    // 第二个单元的节点坐标
    auto points   = SG::Algebra::PointArray (4);
    points[0].m_x = 3.0;
    points[0].m_y = 0.0;
    points[0].m_z = 0.0;
    points[1].m_x = 3.0;
    points[1].m_y = 1.011;
    points[1].m_z = 0.0;
    points[2].m_x = 2.0;
    points[2].m_y = 1.0;
    points[2].m_z = 0.0;
    points[3].m_x = 2.0;
    points[3].m_y = 0.0;
    points[3].m_z = 0.0;
    ShearElementCalculator elementCalc;
    Real                   area;
    Real                   warp;
    Matrixd                Tel (3, 3);

    SG::Algebra::PointArray localPoints2 (4);
    elementCalc.ComputeCoordTrans (points, Tel, localPoints2, area, warp);

    Real pSet[4];
    auto adaptability2 = elementCalc.computeAdaptability (localPoints2, area, G2, property2, nu2, pSet);

    Real adaptabilityTarget2 = 1.778518068655196e-05;
    ASSERT_NEAR (adaptability2, adaptabilityTarget2, SG::DataStructure::Common::TOLERANCE);

    Matrixd Kel2 (12, 12);
    Matrixd Tel23 (2, 3);
    Tel23 (0, 0) = Tel (0, 0);
    Tel23 (0, 1) = Tel (0, 1);
    Tel23 (0, 2) = Tel (0, 2);
    Tel23 (1, 0) = Tel (1, 0);
    Tel23 (1, 1) = Tel (1, 1);
    Tel23 (1, 2) = Tel (1, 2);
    Matrixd T (8, 12);
    T.BlockFill (0, 0, Tel23);
    T.BlockFill (2, 3, Tel23);
    T.BlockFill (4, 6, Tel23);
    T.BlockFill (6, 9, Tel23);
    elementCalc.ComputeStiffness (localPoints2, adaptability2, area, E2, property2, Kel2);
    Kel2 = T.Transpose () * Kel2 * T;

    Matrixd Kel2Target(12, 12, {
        7028.21429000000	,-7028.21429000000	,0	,-7028.21429000000	,-7105.52465000000	,0		,-7028.21429000000	,7028.21429000000	,0		,7028.21429000000	,7105.52465000000	,0,
        -7028.21429000000	,7028.21429000000	,0	,7028.21429000000	,7105.52465000000	,0		,7028.21429000000	,-7028.21429000000	,0		,-7028.21429000000	,-7105.52465000000	,0,
        0	                ,0	                ,0	,0	                ,0	                ,0		,0	                ,0                  ,0		,0	                ,0                	,0,
        -7028.21429000000	,7028.21429000000	,0	,7028.21429000000	,7105.52465000000	,0		,7028.21429000000	,-7028.21429000000	,0		,-7028.21429000000	,-7105.52465000000	,0,
        -7105.52465000000	,7105.52465000000	,0	,7105.52465000000	,7183.68542000000	,0		,7105.52465000000	,-7105.52465000000	,0		,-7105.52465000000	,-7183.68542000000	,0,
        0	                ,0	                ,0	,0	                ,0	                ,0		,0	                ,0                  ,0		,0	                ,0                  ,0,
        -7028.21429000000	,7028.21429000000	,0	,7028.21429000000	,7105.52465000000	,0		,7028.21429000000	,-7028.21429000000	,0		,-7028.21429000000	,-7105.52465000000	,0,
        7028.21429000000	,-7028.21429000000	,0	,-7028.21429000000	,-7105.52465000000	,0		,-7028.21429000000	,7028.21429000000	,0		,7028.21429000000	,7105.52465000000	,0,
        0	                ,0	                ,0	,0	                ,0	                ,0		,0	                ,0                 	,0		,0	                ,0                 	,0,
        7028.21429000000	,-7028.21429000000	,0	,-7028.21429000000	,-7105.52465000000	,0		,-7028.21429000000	,7028.21429000000	,0		,7028.21429000000	,7105.52465000000	,0,
        7105.52465000000	,-7105.52465000000	,0	,-7105.52465000000	,-7183.68542000000	,0		,-7105.52465000000	,7105.52465000000	,0		,7105.52465000000	,7183.68542000000	,0,
        0	                ,0                  ,0  ,0	                ,0                  ,0	    ,0                  ,0                	,0	    ,0                 	,0                  ,0
    });
    ASSERT_TRUE (Kel2 == Kel2Target);
}

// 1,2节点边平行 3，4节点平行边
TEST_F (ShearElementCalculatorTest, ComputeStiffnessS12ParallS34)
{
    auto points   = SG::Algebra::PointArray (4);
    points[0].m_x = 1.0;
    points[0].m_y = -2.0;
    points[0].m_z = 0.0;
    points[1].m_x = 1.11;
    points[1].m_y = -1.0;
    points[1].m_z = 0.0;
    points[2].m_x = 0.0;
    points[2].m_y = -1.0;
    points[2].m_z = 0.0;
    points[3].m_x = 0.0;
    points[3].m_y = -2.0;
    points[3].m_z = 0.0;
    ShearElementCalculator elementCalc;
    Real                   area;
    Real                   warp;
    Matrixd                Tel (3, 3);

    SG::Algebra::PointArray localPoints4 (4);
    elementCalc.ComputeCoordTrans (points, Tel, localPoints4, area, warp);

    Real pSet[4];
    auto adaptability4 = elementCalc.computeAdaptability (localPoints4, area, G2, property2, nu2, pSet);

    Real adaptabilityTarget4 = 1.8775650042624806e-05;
    ASSERT_NEAR (adaptability4, adaptabilityTarget4, SG::DataStructure::Common::TOLERANCE);

    Matrixd Kel4 (12, 12);
    Matrixd Tel23 (2, 3);
    Tel23 (0, 0) = Tel (0, 0);
    Tel23 (0, 1) = Tel (0, 1);
    Tel23 (0, 2) = Tel (0, 2);
    Tel23 (1, 0) = Tel (1, 0);
    Tel23 (1, 1) = Tel (1, 1);
    Tel23 (1, 2) = Tel (1, 2);
    Matrixd T (8, 12);
    T.BlockFill (0, 0, Tel23);
    T.BlockFill (2, 3, Tel23);
    T.BlockFill (4, 6, Tel23);
    T.BlockFill (6, 9, Tel23);
    elementCalc.ComputeStiffness (localPoints4, adaptability4, area, E2, property2, Kel4);
    Kel4 = T.Transpose () * Kel4 * T;

    Matrixd Kel4Target(12, 12, {
        6657.55911000000	,-6657.55911000000	,0		,-7389.89061000000	,-6657.55911000000	,0		,-6657.55911000000	,6657.55911000000	,0		,7389.89061000000	,6657.55911000000	,0,
        -6657.55911000000	,6657.55911000000	,0		,7389.89061000000	,6657.55911000000	,0		,6657.55911000000	,-6657.55911000000	,0		,-7389.89061000000	,-6657.55911000000	,0,
        0	                ,0	                ,0		,0                  ,0                	,0		,0	                ,0	                ,0		,0                	,0	                ,0,
        -7389.89061000000	,7389.89061000000	,0		,8202.77858000000	,7389.89061000000	,0		,7389.89061000000	,-7389.89061000000	,0		,-8202.77858000000	,-7389.89061000000	,0,
        -6657.55911000000	,6657.55911000000	,0		,7389.89061000000	,6657.55911000000	,0		,6657.55911000000	,-6657.55911000000	,0		,-7389.89061000000	,-6657.55911000000	,0,
        0	                ,0	                ,0	  	,0                 	,0                  ,0		,0                  ,0                  ,0		,0                	,0                	,0,
        -6657.55911000000	,6657.55911000000	,0		,7389.89061000000	,6657.55911000000	,0		,6657.55911000000	,-6657.55911000000	,0		,-7389.89061000000	,-6657.55911000000	,0,
        6657.55911000000	,-6657.55911000000	,0		,-7389.89061000000	,-6657.55911000000	,0		,-6657.55911000000	,6657.55911000000	,0		,7389.89061000000	,6657.55911000000	,0,
        0	                ,0                	,0		,0                  ,0                  ,0		,0                  ,0                 	,0		,0                 	,0                	,0,
        7389.89061000000	,-7389.89061000000	,0		,-8202.77858000000	,-7389.89061000000	,0		,-7389.89061000000	,7389.89061000000	,0		,8202.77858000000	,7389.89061000000	,0,
        6657.55911000000	,-6657.55911000000	,0		,-7389.89061000000	,-6657.55911000000	,0		,-6657.55911000000	,6657.55911000000	,0		,7389.89061000000	,6657.55911000000	,0,
        0	                ,0                  ,0		,0                 	,0	                ,0		,0                  ,0                 	,0		,0                  ,0                	,0
    });
    ASSERT_TRUE (Kel4 == Kel4Target);
}

TEST_F (ShearElementCalculatorTest, ComputeStiffnessArbitary)
{
    // 单元的节点坐标
    auto points   = SG::Algebra::PointArray (4);
    points[0].m_x = 3.0;
    points[0].m_y = -2.0;
    points[0].m_z = 0.0;
    points[1].m_x = 3.11;
    points[1].m_y = -0.89;
    points[1].m_z = 0.0;
    points[2].m_x = 2.0;
    points[2].m_y = -1.0;
    points[2].m_z = 0.0;
    points[3].m_x = 2.0;
    points[3].m_y = -2.0;
    points[3].m_z = 0.0;
    ShearElementCalculator elementCalc;
    Real                   area;
    Real                   warp;
    Matrixd                Tel (3, 3);

    SG::Algebra::PointArray localPoints5 (4);
    elementCalc.ComputeCoordTrans (points, Tel, localPoints5, area, warp);
    Real pSet[4];
    auto adaptability5 = elementCalc.computeAdaptability (localPoints5, area, G2, property2, nu2, pSet);

    Real adaptabilityTarget5 = 2.4342817195847756e-05;
    ASSERT_NEAR (adaptability5, adaptabilityTarget5, SG::DataStructure::Common::TOLERANCE);

    Matrixd Kel5 (12, 12);
    Matrixd Tel23 (2, 3);
    Tel23 (0, 0) = Tel (0, 0);
    Tel23 (0, 1) = Tel (0, 1);
    Tel23 (0, 2) = Tel (0, 2);
    Tel23 (1, 0) = Tel (1, 0);
    Tel23 (1, 1) = Tel (1, 1);
    Tel23 (1, 2) = Tel (1, 2);
    Matrixd T (8, 12);
    T.BlockFill (0, 0, Tel23);
    T.BlockFill (2, 3, Tel23);
    T.BlockFill (4, 6, Tel23);
    T.BlockFill (6, 9, Tel23);
    elementCalc.ComputeStiffness (localPoints5, adaptability5, area, E2, property2, Kel5);
    Kel5 = T.Transpose () * Kel5 * T;

    Matrixd Kel5Target(12, 12, {
        6326.81496000000	,-6326.81496000000	,0		,-6953.79662000000	,-6953.79662000000	,0		,-6326.81496000000	,6326.81496000000	,0		,6953.79662000000	,6953.79662000000	,0,
        -6326.81496000000	,6326.81496000000	,0		,6953.79662000000	,6953.79662000000	,0		,6326.81496000000	,-6326.81496000000	,0		,-6953.79662000000	,-6953.79662000000	,0,
        0	                ,0                  ,0		,0                  ,0                 	,0		,0                 	,0	                ,0		,0                 	,0                	,0,
        -6953.79662000000	,6953.79662000000	,0		,7642.91160000000	,7642.91160000000	,0		,6953.79662000000	,-6953.79662000000	,0		,-7642.91160000000	,-7642.91160000000	,0,
        -6953.79662000000	,6953.79662000000	,0		,7642.91160000000	,7642.91160000000	,0		,6953.79662000000	,-6953.79662000000	,0		,-7642.91160000000	,-7642.91160000000	,0,
        0	                ,0                  ,0		,0                 	,0	                ,0		,0                 	,0                  ,0		,0                  ,0                	,0,
        -6326.81496000000	,6326.81496000000	,0		,6953.79662000000	,6953.79662000000	,0		,6326.81496000000	,-6326.81496000000	,0		,-6953.79662000000	,-6953.79662000000	,0,
        6326.81496000000	,-6326.81496000000	,0		,-6953.79662000000	,-6953.79662000000	,0		,-6326.81496000000	,6326.81496000000	,0		,6953.79662000000	,6953.79662000000	,0,
        0	                ,0                	,0		,0                 	,0                 	,0		,0                 	,0	                ,0		,0                  ,0                	,0,
        6953.79662000000	,-6953.79662000000	,0		,-7642.91160000000	,-7642.91160000000	,0		,-6953.79662000000	,6953.79662000000	,0		,7642.91160000000	,7642.91160000000	,0,
        6953.79662000000	,-6953.79662000000	,0		,-7642.91160000000	,-7642.91160000000	,0		,-6953.79662000000	,6953.79662000000	,0		,7642.91160000000	,7642.91160000000	,0,
        0	                ,0	                ,0		,0	                ,0	                ,0		,0	                ,0	                ,0		,0	                ,0	                ,0
        });
    ASSERT_TRUE (Kel5 == Kel5Target);
}

TEST_F (ShearElementCalculatorTest, ComputeMass)
{
    ShearElementCalculator  elementCalc;
    Real                    area;
    Real                    warp;
    Matrixd                 Tel (3, 3);
    SG::Algebra::PointArray localPoints (4);
    elementCalc.ComputeCoordTrans (points, Tel, localPoints, area, warp);
    Real    rho = 7.85;
    Matrixd Mass;
    elementCalc.ComputeMass (localPoints, rho, property, Mass);

    Matrixd MassTarget (
        12,
        12,
        { 1.46062762E-02, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 1.46062762E-02, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 1.46062762E-02, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 1.44615992E-02,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 1.44615992E-02, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 1.44615992E-02, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 1.68986226E-02, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 1.68986226E-02, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          1.68986226E-02, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 1.70432996E-02, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 1.70432996E-02, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
          0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 1.70432996E-02 });

    ASSERT_TRUE (Mass == MassTarget);
}

TEST_F (ShearElementCalculatorTest, ComputeForce)
{
    // Input data
    Real  area = 0.455;
    Real  warp = 0;
    Matrixd Tel23   = Matrixd (2, 3, { -0.24253562503633289, 0.97014250014533177, 0, -0.97014250014533177, -0.24253562503633289, 0 });
    Matrixd eleDisp = Matrixd (12,
                               1,
                               {
                                   0.0,
                                   0.0,
                                   0.0,
                                   -1.1933653795018892E-8,
                                   0.0,
                                   0.0,
                                   0.0,
                                   2.3471272132045944E-8,
                                   2.9224093098158896E-6,
                                   2.0124003538361324E-6,
                                   3.0706144617841325E-7,
                                   5.5933060507344225E-6,
                               });
    IsotropicMaterial material;
    material.m_E  = 210000;
    material.m_nu = 0.3;
    material.m_G  = 80769.234;
    PSHEARData property;
    property.m_thickness = 0.35;
    SG::Algebra::PointArray localPoints (4);
    localPoints[0] = Point{0.0, 0.0, 0.0};
    localPoints[1] = Point{0.82462112512353203, 0.0, 0.0};
    localPoints[2] = Point{0.87312825013079853, 0.60633906259083237, 0.0};
    localPoints[3] = Point{0.26678918753996617, 0.65484618759809887,  0.0};

    // Test data
    SHEARForceData eleForceData;
    SHEARStressData eleStressData;
    ShearElementCalculator elementCalc;
    elementCalc.ComputeElementForceAndStress (localPoints, property, material, area, warp, eleDisp, Tel23, _OUT eleForceData, _OUT eleStressData);

    // Target data
    SHEARForceData eleForceDataTarget;
    eleForceDataTarget.m_F41 = 0.011832801338066627;
    eleForceDataTarget.m_F21 = 0.011315425242898523;
    eleForceDataTarget.m_F12 = -0.011315425242898523;
    eleForceDataTarget.m_F32 = -0.012181739866474962;
    eleForceDataTarget.m_F23 = 0.012181739866474962;
    eleForceDataTarget.m_F43 = 0.014855780324969472;
    eleForceDataTarget.m_F34 = -0.014855780324969472;
    eleForceDataTarget.m_F14 = -0.011832801338066627;
    eleForceDataTarget.m_KF1 = 0;
    eleForceDataTarget.m_S12 = 0.027443937338382927;
    eleForceDataTarget.m_KF2 = 0;
    eleForceDataTarget.m_S23 = 0.040053313953315635;
    eleForceDataTarget.m_KF3 = 0;
    eleForceDataTarget.m_S34 = 0.048845504821123001;
    eleForceDataTarget.m_KF4 = 0;
    eleForceDataTarget.m_S41 = 0.03346821626632066;

    // Check
    ASSERT_TRUE (eleForceData == eleForceDataTarget);
}

TEST (ShearElementWarpTest, ComputeForceandStress)
{
    Real area = 86.314639450129334;
    Real warp = 0.63727252791867739;
    Matrixd Tel23   = Matrixd (2, 3, {0.98892809007285765,
                                      0.12928618137949993,
                                      -0.072844464231382827,
                                      -0.12230396088377024,
                                      0.98809631120024732,
                                      0.093313562490162996});
    Matrixd eleDisp = Matrixd (12,
        1,
        {
            -0.21125896562490915,
            0.20706933515931336,
            0.016779095361052928,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.14505064148208743,
            0.46358615082972904,
            -1.0799042116033786,
        });
    IsotropicMaterial material;
    material.m_E  = 210000;
    material.m_nu = 0.29999999999999999;
    material.m_G  = 80769.230769230766;
    PSHEARData property;
    property.m_thickness = 1;
    SG::Algebra::PointArray localPoints (4);
    localPoints[0].m_x = 0.0;
    localPoints[0].m_y = 0.0;
    localPoints[0].m_z = 0.0;
    localPoints[1].m_x = 8.6868850695229138;
    localPoints[1].m_y = 0.0;
    localPoints[1].m_z = 0.0;
    localPoints[2].m_x = 9.8826092050865437;
    localPoints[2].m_y = 7.7147775043799403;
    localPoints[2].m_z = 0;
    localPoints[3].m_x = -0.76015887841607188;
    localPoints[3].m_y = 10.093228893081264;
    localPoints[3].m_z = 0;
    SHEARForceData eleForceData;
    SHEARStressData eleStressData;
    ShearElementCalculator elementCalc;
    elementCalc.ComputeElementForceAndStress (localPoints, property, material, area, warp, eleDisp, Tel23, _OUT eleForceData, _OUT eleStressData);

    SHEARForceData eleForceDataTarget;
    eleForceDataTarget.m_F41 = -1313.9234681459714;
    eleForceDataTarget.m_F21 = -1775.8027802095976;
    eleForceDataTarget.m_F12 = 1775.8027802095976;
    eleForceDataTarget.m_F32 = 1646.1146927094001;
    eleForceDataTarget.m_F23 = -1646.1146927094001;
    eleForceDataTarget.m_F43 = -1461.4907205162738;
    eleForceDataTarget.m_F34 = 1461.4907205162738;
    eleForceDataTarget.m_F14 = 1313.9234681459714;
    eleForceDataTarget.m_KF1 = -212.99846092707503;
    eleForceDataTarget.m_S12 = -408.84684579051884;
    eleForceDataTarget.m_KF2 = -264.6449325372617;
    eleForceDataTarget.m_S23 = -421.70812384417781;
    eleForceDataTarget.m_KF3 = -219.77657770739944;
    eleForceDataTarget.m_S34 = -268.03313465691116;
    eleForceDataTarget.m_KF4 = -168.13010609721272;
    eleForceDataTarget.m_S41 = -259.62214226620728;

    SHEARStressData eleStressDataTarget;
    eleStressDataTarget.m_TAVG = -348.97085891374167;
    eleStressDataTarget.m_TMAX = 509.65082530884092;

    ASSERT_TRUE(eleForceData == eleForceDataTarget);
    ASSERT_TRUE(eleStressData == eleStressDataTarget);
}